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Abstract-In the present paper, the Nonlinear Harmonic (NLH) 
method implemented in Fine/Turbo code has been evaluated 
for unsteady rotor-rotor interaction analysis. Strengths and 
weaknesses of the NLH method are highlighted through direct 
comparison with experimental data and time-domain unsteady 
simulation for the NWPU counter-rotating research 
compressor (CRRC) at the design and off-design conditions. 
The rotor-rotor interactions were studied in detail to advance 
the understanding of the flow physics involved. Based on the 
unsteady simulation, deterministic stresses were computed and 
compared to the NLH method to analyze distinct sources of 
unsteadiness and their effect on the time-averaged flow field. 
For the present case, the NLH method is about one order 
magnitude faster than the direct time-domain method. 
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I. INTRODUCTION 

The current generation of high performance fans and 
compressors are being designed with the goal of 
maximizing blade loading and minimizing blade row 
spacing in order to increase performance and reduce 
compressor length and thus weight. Counter-rotation, as one 
of the promising technologies in lowering engine weight, 
has become more and more attractive for substantial 
improvement in aerodynamic performance. By counter 
rotation the pressure ratio can be increased substantially 
without needing high circumferential speed in the individual 
blade rows because the sumof the circumferential velocities 
of two adjacent blade rows is available for the work input. 
The stator could be cancelled from the compressor and the 
torque of the two shafts in a counter-rotating system is less 
than half of that in a single shaft system The result is a 
potential for higher pressure ratios with fewer blade rows, 
hence either shorter or lighter compressors. Applications for 
counter-rotating technology do have the potential to be an 
alternative propulsor for future very high bypass engines 
and in high-supersonic cruise aircraft [2 3] . 

E. NONLINEAR HARMONIC METHOD 

The time-domain method to solve unsteady flow fields 
in multistage turbomachines can provide insights to help our 
understanding of complex unsteady flow phenomena. 
However, due to unequal blade counts in neighboring 
bladerows, an unsteady calculation has to be carried out in 
multiple passages/whole annulus owing to the required 
direct periodic condition. Moreover, it is recognized that the 
time-domain unsteady simulation is difficult to extract 
useful information and prohibitively expensive in 



computational time, which, in an engineering context, can 
constitute a major limitation. Therefore it is not feasible to 
be used in a routine design and explains why steady 
simulation is still widely used at industrial framework. 

The development of efficient unsteady aerodynamic 
modeling methodologies has been an active research area. In 
order to overcome that problem, the sophisticated frequency 
domain approach has been extended into the multistage 
turbomachiny to predict blade row interactions. He[4] 
proposed a nonlinear harmonic (NLH) method that revealed 
superior to the time -domain unsteady simulation. 

The main idea in NLH method is that an unsteady flow 
variable can be decomposed into a time -averaged part and a 
sum of periodic perturbations, which in turn can be 
decomposed into N harmonics. 

U(rj) = U(r) + YF'( 7 > t ) (1) 

Where U(f) is the time-averaged part, U\r,t) is the 
unsteady part induced by the unsteady disturbance. The 
unsteady disturbances are caused by blade passing 
frequency of the adjacent rotating blade row. These 
periodically perturbations are transformed into the 
frequency domain by a Fourier transform, which represents 
the oscillating influence of the perturbations caused by the 
adjacent rows. The accuracy of the solution is determined 
by the number of harmonics and, especially for multi-stage 
environment, the number of perturbations or blade passing 
frequencies. It is highly desired that this algorithm solve for 
the periodic solution directly, avoiding the expensive 
transients typical of time -domain schemes and consume 
most of the CPU resources. 

The so-called deterministic stress terms due to 
nonlinearity of the momentum and energy equations were 
solved simultaneously with the harmonic perturbation 
equations in a strongly coupled approach. These terms are 
similar to the Reynolds stress terms, which can account for 
the nonlinear effects of unsteadiness on the time-averaged 
flow, making this approach similar to the average-passage 
method of Adamczyk [5] . The important advantage of the 
method is that the extra deterministic stress is calculated 
directly from products of the real parts and imaginary parts 
of the complex amplitudes [6] . 

HI. FOCUS OF RESEARCH EFFORT 
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The NLH method has been tested and validated against 
several turbomachinery cases, with the result that the time 
savings seemed to be significant over the time -do main 
unsteady method [6 7] . The first priority of present study is 
directed toward demonstrating those claims in a counter- 
rotating research compressor, as the NLH method is 
relatively new with limited validation. Apart from this, there 
are several issues of interest. First, the performance maps of 
NLH method results are verified by comparison with 
experimental data. Second, considering that the validation of 
a 3-D unsteady flow structure is much more difficult than its 
steady counterpart because of the extreme difficulty to 
obtain detailed unsteady data to provide quantitative 
measurement and visualization of flow structure in CRRC. 
Therefore, to demonstrate the accuracy and computational 
efficiency of the NLH method, it is also compared to a time- 
domain unsteady simulation to advance our present 
understanding of the complicated interactions between 
rotors, which is deemed to be a real challenge to the NLH 
method. 

IV. EXPEREVI ENT A L FA CILIT Y 

The numerical test rig in this work is a counter-rotating 
research compressor (CRRC) at the NWPU National 
Defense Aerodynamics Laboratory of Airfoil and Cascade, 
which is shown in Fig. 1. Over the past three years we have 
done much research on the CRRC [8 9] . The CRRC consists 
of a variable inlet guide vane (with 20 airfoils), the rotor-I 
(with 19 airfoils), the rotor-II (with 20 airfoils) and the 
outlet guide vane (with 32 airfoils). The design speed ratio 
of the CRRC was chosen to be 1 and the wheel speeds of the 
two counter-rotating rotors were 8000 and -8000 rpm 
respectively. The stage design total pressure ratio is 1.22 at 
a corrected mass flow rate of 6.4Kg/s. Both rotors 
configurations were modeled with the same tip clearance of 
0.5mm. More details about the aerodynamic design and 
experimental measurement are given by Liu [10] and Chen [11] . 




Fig. 1 ViewoftheCAAC 



V. DESCRIPTION OF THE NUMERICAL MODEL 

The numerical simulation were undertaken with the 
commercial CFD code Fine™/Turbo of NUMECA, which 
has been extensively used for turbomachinery industry. Its 



continuous development over the years has extended its 
versatility to a number of aero -engines design applications. 

The geometry of the CRRC has been partitioned using 
pre-processor Autogrid5. A structured O4H mesh was 
generated for each blade passage, which is convenient for 
parallel computation. An O-grid is generated around the 
blade for a good control of the stretching in the direction 
normal to the wall. A butterfly grid topology and 33 grid 
points is applied in the tip clearance region of both rotors to 
model leakage effects accurately. The mesh is refined 
towards the solid boundaries to meet the requirements of the 
low-Reynolds turbulence model of y+ < 5. Since the NLH 
method works in the frequency domain, only one blade 
passage is needed. The number of grids in each of the rotor 
passages is listed in Table I. 

Past simulations of Vilmin [12] indicated that the number 
of grid points in the circumferential direction has to follow 
some strict guidelines. He suggested that there be at least 30 
points in the circumferential direction for each harmonic. 
Furthermore, the number of grid points has to be weighted 
according to the number of blades on both sides of the 
rotor/rotor interface. Following his suggestions, the number 
of grid points in each of the blade passages is listed in Table 
1. This grid density was deemed adequate for the purpose of 
this study. 

Common grid has been subsequently used for the time- 
domain unsteady simulation to ensure consistency. While 
the original blade row number ratio is 20:19:20:32, for the 
time-domain unsteady simulation, it was necessary to 
change the blade row number ratio to 20:20:20:30. This 
simplification afforded a 1/1/1/2 model, which is a 
considerable reduction in computational resources. A 
modicum of scaling was used for rotor I and OGV to 

preserve the same solidities so that the blade loadings were 
n 31 

approximately constant 

Prior to the NLH calculation, a single passage steady 
calculation with blade row interface mixing plane is carried 
out. For the steady flow, the governing equations are 
discretized in space using a cell-centered finite volume 
scheme together with the blend second-order and fourth- 
order artificial dissipation to damp numerical oscillations. 
Temporal integration of the governing equations is carried 
out using the explicit four-step Runge-Kutta scheme. The 
calculation is accelerated using local time -stepping and 
multi-grid techniques. The Reynolds stresses are closed by 
the S-A turbulence model because of its excellent stability. 

At the inlet plane the measured total pressure and 
temperature quantities are imposed while the measured 
average static pressure was imposed at the exit plane. Non- 
slip and adiabatic conditions are adopted for all the solid 
walls. At blade row interface, unsteady disturbances for 
each blade row are obtained by a Fourier transform of 
pitch wise spatial nonuniformity of an adjacent blade row. 
The ID non-reflecting boundary conditions are applied in 
order to prevent spurious reflections of outgoing waves back 
into the computational domain. The phase-shifted boundary 
conditions are used on the upper and lower periodical 
boundaries to solve the unsteady perturbation equations in a 
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single passage. It is essentially equivalent to the 
chorochronic conditions demonstrated by Gerolymos [14] . In 
the present simulation, the downstream running wake and 
upstream running potential wave is resolved by the first 
three harmonics. 

The main drawback of the NLH method is that it 
requires much more memory compared to a typical time- 
domain unsteady simulation [12] . Considering this and in 
order to reduce computation time, we choose a cluster 
system to do the coarse-grained parallel simulation, which 
consists of 14 IBM BladeCenter HS21 nodes and runs in a 
Linux environment. Each node has 8GB memory size and 2- 
way 2 GHz Intel Xeon E5405 Quad core processors. 
Parallel Virtual Machine (PVM) and Open Message Passing 
Interface (OM PI) communication protocols are used for 
Interprocess or communication. 

The time -do main unsteady method is based on the dual- 
time stepping approach of Jameson [19] , where every 
instantaneous flow solution is considered as a "pseudo- 
steady" state problem 40 physical time steps with 40 inner 
pseudo-time steps per Rotor 2 b lade passing are used. 

The entire operating range of the CRRC, from choke to 
near-stall conditions, is simulated with nine different mass 
flow rate values. Efficiency, mass flow, pressure ratio and 
residuals were monitored during the computations for 
assuring that the solution was converged. In terms of 
computational efficiency, the NLH method starting from a 
converged steady solution could be obtained in about 20-30 
hours depending on the simulation point. This represents a 
big savings in computational time while the time-domain 
unsteady calculation needs about 20 hours for one period 
and at least 10 periods to obtain a periodic solution. 
Therefore, for this case, the NLH method is about one order 
magnitude faster than the direct time-domain method. 

TABLE I NUMBER OF GRIDS FOR THE BLADE PASSA GE 
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VI. COMPARISON WITH EXPERIMENTAL DATA 

Figure 2 shows the experimentally measured 
performance maps along with the numerical results obtained 
using the present NLH method and time-average of 
unsteady simulation. The experimental values are area 
averaged quantities while the numerical simulations are 
mass averaged quantities. It can be seen that the total 
pressure ratio of both simulations is in good agreement with 
measured results, while isentropic efficiency is a little 
deviation from simulations. As analyzed by Chen, the 
measured efficiency have a measurement uncertainty which 
is dominated by temperature rise uncertainty [8] . It should be 
noted that the mass flow rate resulting from NLH simulation 



is larger than that from time-averaged unsteady simulation 
for the whole operating line. In design point the efficiency 
predicted by the NLH method is only 0.5% higher than that 
predicted by the time-averaged unsteady simulation. When 
the operating condition moves away from the design point, 
the discrepancy between them grows larger and is about 3% 
for the near stall condition. The discrepancy might stem 
from separated flow fields at near-stall condition, which 
exhibits strong nonlinear unsteady interactions and cannot 
be resolved with the NLH method with only three orders. As 
analyzed by He [16], conventional harmonic method using a 
pseudo -time- marching technique often exhibits solution 
divergence behavior for separated flow fields while it is 
often unavoidable for highly loaded design at near-stall 
conditions. In total, the satisfying results can be achieved 
with the NLH method in a fraction of time in contrast to a 
time-domain unsteady simulation although small 
discrepancies exist, also the NLH results are very accurately 
reproduced compared with the test data. 

Figures 3 and 4 show profiles of the radial distribution 
of the mass -averaged absolute total pressure and isentropic 
efficiency at the exit plane of the CRRC for NLH method 
and unsteady simulation. The agreement between the two 
sets of results is also very good. Thus the NLH method 
provides a correct representation of the time-averaged flow 
field to the CRRC. 




5 5.5 S 65 7 7 5 6 
Mass Flow (Kg/S) 

Fig.2 Comparison of experiment, NLH method and time-averaged unsteady 
characteristic 




20 40 80 1CO 

Span % 



Fig. 3 Span wise distribution of isentropic efficiency 
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Fig. 4 Span wise distribution of total pressure ratio 

VII. INFLUENCE OF ROTOR-ROTOR 
INTERACTION 

Comparisons of the NLH method and time-domain 
unsteady method of CRRC with steady simulation are 
carried out at the aero design point. Figure 5 shows the 
comparison of one-dimensional distributions of entropy 
along axial direction in the rotor I and rotor II passage 
respectively, which are area averaged for the whole S3 
plane. The focus here is placed on the rotor/rotor interface. 
In steady simulation, only the averaged flux can be 
transferred through the interface, and the non -uniform flux 
is completely mixed out at the rotor I outlet boundary and 
subsequently transferred to the rotor II frame to provide 
uniform and steady inflow boundary conditions. 
Considerable loss is generated across the interface due to the 
mixing of wake blockage. The mixing loss at interface is 
removed by adding the deterministic stresses in the NLH 
method. It shows that there is little difference between the 
NLH and unsteady simulation in most of the axial positions. 
With three orders of perturbation, the NLH method can 
recover nearly 80% of the mixing loss. Continuity of the 
flow variables across the interface is exactly satisfied for an 
infinite number of time harmonics. Since the NLH method 
uses a limited number of harmonics, this continuity cannot 
be rigorously reproduced. But the observed continuity in the 
numerical results increases with a higher number of 
harmonics . 
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Fig. 5 Comparisons of 1-D entropy distribution along axial direction 



Figure 6 with the inset shows the harmonic static 
pressure amplitude on rotor-II blade surface at 50% span of 
aero design point with first harmonic, second harmonic, and 
third harmonic. The first two orders of perturbations 
dominate the pressure amplitude, which is consistent with 
the previous analysis that the first three orders of unsteady 
perturbations can recover about 80% of the mixing loss. 
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Fig. 6 Harmonic pressure amplitude on rotor-II blade surface at 50% span 

Figure 7 shows the reconstructed instantaneous entropy 
contours at 50% span for rotor I, rotor II and OGV, as 
derived from the NLH method. The time increments are 
equally spaced over a cycle of the rotor-rotor period. The 
efforts in terms of mesh generation and the simulation 
configuration have allowed for a good resolution of the 
important flow phenomena. As can be seen in the figure, the 
flow field exhibits strong unsteady phenomena that are 
periodic in nature. 
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Fig. 7 Instantaneous entropy 

The rotor I and rotor II wakes are clearly visible as 
regions of high entropy. At the given point of time, for 
example, t=l/20T, The wakes from rotor I are chopped at 
the rotor II blade leading edge and divided into two 
segments. The chopped wake segments propagate 
independently of each other within the adjacent rotor blade 
channels with the free stream velocity. Depending on the 
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solidity and speed ratio, several chopped rotor I wakes can 
be found within a single rotor II blade passage. For the 
rotor-rotor blade count and design wheel speeds of the 
CRRC, three to four rotor I wake segments appear within a 
single rotor II blade passage, depending on the point of time. 
The propagation velocity of the upstream wake segments is 
determined by the velocity distribution within the blade 
passages. For this reason, in the front part of the blades, the 
wake segments propagate faster near the suction surface 
than near the pressure surface. This is more evident in the 
OGV pas sages. 

Once the rotor I wake are chopped by the rotor II and 
then propagate downstream inside, the low velocity inside 
the wake causes the low momentum wake to transport and 
collect near the pressure surface forming so called 
"negative jet" when the wake and entropy is transported 
downstream This results in a local broadening of the wake 
segment near the pressure surface and a thinning near the 
suction surface. The characteristics of "negative jet' is 
similar to a typical axial-flow compressor [21] . The 
magnitude of wake -induced unsteadiness depends on the 
spacing between the two rotors. With three orders of 
perturbation, the NLH method successfully captures the 
physical phenomenon of "wake stretching" inside the rotor 
II and OGV passage. It compares very well with the 
instantaneous entropy contours produced from the time- 
domain unsteady simulation (not shown). 

Figures 8 and 9 show the rotor I and rotor II surface 
static pressure distributions at 50% span of aero design point. 
Time -averaged and reconstructed instantaneous values 
based on the phase-shifted periodicity are shown on the 
same plots. The time-averaged loading of rotor I is almost 
the same as the instantaneous loading except in the trailing 
edge. In design condition, the relative Mach number of the 
full span of the rotor II are subsonic, so there is no bow 
shock interaction with the upstream rotor wake, which 
means that pressure fluctuations on the rotor I are dominated 
by the rotor II potential fields as they propagate upstream 
the rotor I. 
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Fig. 8 Pressure fluctuations on rotor I surface 
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Fig. 9 Pressure fluctuations on rotor II surface 
These fluctuations appear weak in comparison with the 
fluctuations of the loading on the rotor II. For the rotor II, 
larger fluctuations are observed, which is due to the impact 
of the wake of rotor I propagated downstream and the 
potential fields of OGV propagated downstream inside the 
rotor II passage. The wake -induced unsteadiness and 
circumferential variation of flow conditions in the 
rotor/rotor interface causes period fluctuations of inlet flow 
conditions for the rotor II. Overall, the rotor-rotor 
interaction can affect a range of 85% rotor I axial chord 
from trailing edge performance, while whole axial chord 
length on rotor II performance, which indicates that the 
interaction between rotor I and rotor II is significant. 

VIE. DETERMINISTIC STRESSES STUDY 

The deterministic stresses appear in the NLH method 
through time averaging of the periodic unsteadiness due to 
the unsteady interaction of the neighboring blades. The 
deterministic stress has been post-processed from NLH 
simulation and compared to the time-domain unsteady 
results. An averaging procedure based on the method of 
Hall [18] has been developed to determine the deterministic 
stresses from the time -averaged flow field of the unsteady 
simulation. A total of 40 snapshots in one blade passing 

period were processed. In these figures, only (pu'^ju' are 
presented due to the paper length limitation 

Figure 10(a) shows the full deterministic stresses 
(pv'^u' (the sum of three orders) for NLH method at 

50% span for rotor I. It can be seen that at the inlet the 
deterministic stresses levels are very low, primarily because 
the axial spacing between IGV and rotor I is relatively large 
so the wake of IGV decay rapidly and totally mixed out at 
the interface. As the flow pass through the blade passage, 
higher levels begin to appear near the trailing edge, as a 
result of potential effects from the interaction with the 
downstream rotor. In circumferential direction, the 
deterministic stresses are not uniform and the gradients are 
large in some regions. This is contrast to the average 
passage equation proposed by Adamczyk [5] where the 
deterministic stresses are uniform across the circumferential 
direction because it is rather hard and time-consuming to 
take account for the circumferential distributions in the 
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blade passages. In Figure 10 (b) Fairly good agreement is 
obtained based on unsteady simulation compared to NLH 
method. 

Figures 10(c) present the results for NLH method at 
50% span for rotor II. It can be seen that the deterministic 
stresses in the rotor I and rotor II have the same 
magnitudes .The location of high concentrations of 

(pu'^L)' is observed on the rotor/rotor interface. The 

deterministic stresses reduce in strength as the flow moves 
downstream as the unsteadiness due to the Rotor I wake 
dissipates. As the flow passes through the blade passage, 
higher levels begin to appear in the rotor wake, presumably 
as a result of potential effects from the interaction with the 
downstream OGV. The unsteady simulation results present 
also high values upstream of the rotor II leading edge and 
the transport and development of the deterministic stress in 
the flow field, which shows a good match between the two 
methods. Also the comparisons have been achieved for all 
of the other deterministic stresses (not shown). 

Figures 11 present the deterministic stresses for 99% 
span. The patterns of the predicted deterministic stresses 
near the tip regions were considerably different from the 
mid-span where the tip leakage vortex is instead the 
dominant source of unsteadiness, especially in rotor II. One 
of the striking features is that a significant source of high 
deterministic stresses begins roughly from the blade leading 
edge following the approximate trajectory of the tip 
clearance flow. The NLH method manages to reproduce the 
positions of the local extrema of the deterministic stresses in 
rotor I. However, in the leading edge of rotor II there is a 
noticeable difference in the quantitative results, partly 
because of the strong upstream wakes and rotor II tip 
clearance flow interaction, which has ever higher blade 
passing frequency and needs more orders of harmonics to 
simulate accurately. 
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Fig. 10 comparison of deterministic stress 
( po'^u' distribution at 50% span 
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Fig. 11 comparison of deterministic stress 

( pu'^ju' distribution at 99% span 
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IX. CONCLUDING REMARKS 

To investigate the capability of the NLH method on 
predicting unsteady flows induced by rotor-rotor 
interactions, the unsteady flows in a counter-rotating 
research compressor (CRRC) have been simulated. Starting 
from a steady solution, a converged solution with 3 
harmonics could be obtained in about 20 hours on IBM 
BladeCenter HS21 using OMPI. The NLH method delivered 
accurate results at design points compared to time-domain 
unsteady results. 

Analysis of the results provided some insight into the 
flow mechanism of rotor-rotor interaction, the following 
conclusions were drawn: 

1. The NLH method takes at least one order of 
magnitude faster than typical time-domain unsteady 
simulations in terms of computational efficiency, which 
seems to be very promising for routine unsteady simulations 
in the design of modern multistage turbomachinery. 

2. Comparative studies show that the present NLH 
method can capture the time -averaged results in good 
agreement with the time -domain unsteady method and 
experimental results, especially for the aero design point. 

3. With the three orders of wake and potential 
perturbations included, 90% mixing loss could be recovered 
at the rotor/rotor interface. 

4. The NLH method has proven to resolve and sustain 
rotor-rotor interaction phenomena quite well. In the present 
test case, the primary contributor to the deterministic 
stresses at mid-span is the interaction of a blade with the 
upstream rotor wakes while the tip leakage flow is the 
dominant source of unsteadiness and deterministic stresses 
in the tip region. 

The achieved results are a valuable enlargement of a 
data base in order to design counter rotating compressor 
with higher performance. 
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